AN OPTIMAL ORDER ERROR ANALYSIS OF THE ONE-DIMENSIONAL 
QUASICONTINUUM APPROXIMATION 



MATTHEW DOBSON AND MITCHELL LUSKIN 

Abstract. We derive a model problem for quasicontinuum approximations that allows a simple, 
yet insightful, analysis of the optimal-order convergence rate in the continuum limit for both the 
energy-based quasicontinuum approximation and the quasi-nonlocal quasicontinuum approxima- 
tion. The optimal-order error estimates for the quasi-nonlocal quasicontinuum approximation are 
given for all strains up to the continuum limit strain for fracture. The analysis is based on an explicit 
treatment of the coupling error at the atomistic to continuum interface, combined with an analysis 
of the error due to atomistic and continuum schemes using the stability of the quasicontinuum 
approximation. 



1. Introduction 

The quasicontinuum method (QC) denotes a class of related approximations of fully atomistic 
models for crystalline solids that reduces the degrees of freedom necessary to compute a deformation 
to a desired accuracy [6, 8, 9, 12, 13, 15-19, 21, 22, 25, 27, 28]. The derivation of a quasicontinuum 
method first removes atomistic degrees of freedom by using a piecewise linear approximation of the 
atom deformations with respect to a possibly much smaller number of representative atoms. Since 
atoms interact significantly with atoms beyond their nearest neighbors, a further approximation is 
required to obtain a computationally feasible method. 

In this paper, we will analyze two QC variants that approximate the total atomistic energy 
by using a continuum approximation in a portion of the material called the continuum region. 
The deformation gradient is assumed to be slowly varying in the continuum region, making the 
continuum approximation accurate. The more computationally intensive atomistic model is used for 
the remainder of the computational domain, which is called the atomistic region. In this region, all 
of the atoms are representative atoms, so that there is no restriction of the types of deformations in 
the atomistic region. To maintain accuracy, the atomistic region must contain all regions of highly 
varying deformation, such as material defects. Adaptive methods that determine what portion of 
the domain should be assigned to the atomistic region in order to acheive the required accuracy 
have been considered in [1-3,20-22,24]. Other approaches to atomistic to continuum coupling have 
been developed and analyzed in [4,23], for example. 

In Section 2, we derive a model problem for QC approximations and describe the energy-based 
quasicontinuum (QCE) approximation and the quasi-nonlocal quasicontinuum (QNL) approxima- 
tion. These two approximations use the same continuum approximation, but differ in how they 



Date: January 29, 2009. 

2000 Mathematics Subject Classification. 65Z05,70C20. 

Key words and phrases, quasicontinuum, error analysis, atomistic to continuum. 

This work was supported in part by DMS-0757355, DMS-0811039, the Institute for Mathematics and Its Applica- 
tions, the University of Minnesota Supercomputing Institute, and the University of Minnesota Doctoral Dissertation 
Fellowship. This work is also based on work supported by the Department of Energy under Award Number DE- 
FG02-05ER25706. 



AN OPTIMAL ORDER ERROR ANALYSIS OF THE QUASICONTINUUM APPROXIMATION 



2 



couple the atomistic and continuum regions. We also give stability estimates for the two QC 
variants. 

We have derived our model quasicontinuum energy from a general quasicontinuum energy by 
expanding each interaction to second order to be able to present a simple, but illuminating, analysis. 
The model differs from a standard quadratic approximation by keeping certain first-order terms. 
These are a source of leading-order coupling error, and reflect the behavior in the non-linear case. 
We have also chosen to analyze the model problem for boundary conditions given by restricting to 
periodic displacements to maintain the simplicity of the analysis. 

The goal of this paper is to give an error analysis as the number of atoms per interval increases 
(the continuum limit). The residual at atoms in the coupling interface is lower order (order 0(l/h) 
and 0(1) for QCE and QNL, respectively) than the residual in either the atomistic or continuum 
region (0(h?) in all cases). However, the corresponding error depends primarily on the sum of the 
residual at the atoms in the atomistic to continuum coupling interface, and this sum has the higher 
order 0(h) due to the cancellation of the lowest order terms when the residual is summed across 
the interface. 

In Section 3, we split the residual for the QCE approximation into the part due to the continuum 
approximation and the part due to coupling the atomistic and continuum regions. The stability of 
the QCE approximation and the 0(h 2 ) estimate for the corresponding residual combine to give an 
optimal order bound on the error due to the continuum approximation. We then derive an explicit 
representation of the coupling error, and we observe that this error is small and decays away from 
the interface since the coupling residual is oscillatory, as described in the preceding paragraph. The 
coupling error is the leading order term in the total error, dominating the continuum approximation 
error. 

We show that the displacement converges at the rate 0(h) in the discrete l°° norm and the rate 
O(h^P) in the 

w i,p norms where h is the interatomic spacing. Our analysis extends the results of 
E, Ming, and Yang [14] that show that the error is 0(1) in the w 1,co norm for the QCE method 
applied to a problem with harmonic interactions and Dirichlet boundary conditions. 

In Section 4, we present the same analysis in the QNL case. Here we show that the improved 
order of accuracy in the coupling interface serves to nearly balance the order of the error due to the 
continuum approximation, and we are consequently able to give higher order optimal error estimates 
for the QNL approximation than for the QCE approximation. We show that the displacement now 
converges at the rate 0(h 2 ) in the discrete l°° norm and the rate 0(h 1+1 / p ) in the w 1,p norms where 
h is the interatomic spacing. E, Ming, and Yang [14] have obtained 0(h) esimtates in the w 1,co 
norm for the Lennard-Jones potential and strains bounded away from the continuum limit strain 
for fracture. We have obtained optimal order error estimates for the discrete l°° and w l,p norms 
for all strains up to the continuum limit strain for fracture. 

This paper extends our analysis of the effect of atomistic to continuum model coupling on the 
total error in the energy-based quasicontinuum approximation [10] to include external forcing. We 
also now include an analysis of the quasi-nonlocal approximation. 



2. One-Dimensional, Linear Quasicontinuum Approximation 

We consider the periodic displacement from a one-dimensional reference lattice with spacing 
h = 1/N, and we denote the positions of the atoms in the reference lattice by 



Xj := jh, 



-co < j < CO. 
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We will derive and analyze the linearization about a uniform deformation gradient F given by the 
deformation 

Vj :=jFh = ja, -oo < j < oo, (2.1) 

which is a lattice with spacing a := Fh. We will then consider perturbations Uj of the lattice yj 
which are 2N periodic in j, that is, we will consider deformations yj where 

Vj ■= yf + uj, -oo < j < oo, 

for 

Uj+2N = Uj, — oo < J < oo. (2.2) 

We will often describe the perturbations Uj satisfying (2.2) as displacements (which they are if yj 
is considered the reference lattice). We thus have that the deformation satisfies 

Uj+2N = Vj + 2F, -oo<j<oo. (2.3) 

We note that neither the reference lattice spacing h nor the uniform lattice spacing a need be the 
equilibrium lattice constant or the well of the interatomic potential. 



2.1. Notation. Before introducing the models, we fix the following notation. We define the back- 
ward differentiation operator, Du, on periodic displacements by 

, Uj — Uj-l 

(Du)j := — — j-* — for — oo < j < oo. 

Then (Du)j is also 2N periodic in j. We will use the short-hand (Du)j = Duj. 
For periodic displacements, u, we define the discrete norms 

N 




\u\\gp := I h ^2 \ u j\ P I > 1 < p < oo, 



-N+l 

HulLoo := max \v,j\, 
l h ~N+l<j<N J 

By including the whole period, we ensure that they are all norms (in particular, ||u|L P = implies 

h 

u = 0). We sum over a single period to make the norms finite. We will also consider periodic 
functions u(x) : E — ► R satisfying 

u(x + 2) = u(x) for (2.4) 

We define corresponding continuous norms 

•l \ i/p 



\u\ 



LP 



( f Y 

Ij \u(x)\ p dx j , 1 < p < oo, 



||n|| Loo := ess sup |«(x)|. 
xe(-i,i) 

We let v! denote the weak derivative of the periodic function u. We note that if Hu'H^p < oo, then 
u(x) is continuous for all x in E and u(— 1) = u(l). We will similarly denote higher order weak 
derivatives of the periodic function u as u" , u'", and u^. 
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2.2. Atomistic Model. We first consider the total energy per period 

£ tot ' h (y) :=£ a ' h {y)-Hy), 
for deformations y satisfying (2.3) where the total atomistic energy per period is 



j=-N+l 



h 



+ 



h 



(2.5) 



(2.6) 



for a two-body interatomic potential <f> (assumptions on the potential are given in Section 2.5), and 
where the total external potential energy per period is 



N 
=-N+l 



(2.7) 



0. 



for periodic dead loads f such that /j+2jv = fj and Ylf=-N+i fj 

We have scaled the atomistic energy per bond in (2.6) by h<j>(r/h) and the external force per 
atom by hfo in (2.7). This scaling permits a continuum limit as h — > 0. If y, f € C°°(R) satisfy 
y{x + 2) = y(x) + 2F, y'(x) > 0, f(x + 2) = f(x), and f(x) dx = 0; if 4>{r) is locally Lipschitz 
for r £ (0, oo); and if we set yj = y{xj) and fj = f(xj), then the energy per period (2.5) converges 
to [5] 



J 4>{y'( x )) - f( x )y( x ) 



dx 



(2.8) 



as N —> oo (which implies h — > 0), where <p{r) = <p(r) + <p(2r). In the following, we linearize the 
atomistic model which leads to a corresponding linearized continuum model. This paper analyzes 
the convergence of two quasicontinuum approximations to the minimizer of the linearized continuum 
model's total energy. 



2.3. Linearized Atomistic Model. We will henceforth consider the linearized version of the 
above energies while reusing the notation £ a,h and £ tot ' h . The total atomistic energy (2.6) be- 
comes [10] 

N 

£ a > h (u) := h 4 

j=-N+l _ 
+<t>2F 

for displacements, u, satisfying the periodic boundary conditions (2.2). Here 4>' F := (f)'(F),(f>p := 
<fi"(F), (j)' 2F := 4>'(2F), 4>2 F := 4>"{2F), where <p is the interatomic potential in (2.6). We have 
removed the additive constant 2(j){F) + 2<p(2F) from the quadratic expansion of the energy, and 
we will remove the additive constant —h^- = _ N+l fjyJ from .F(u) when computing the external 
potential of the displacement u. We note that the first order terms in (2.9) sum to zero by the 
periodic boundary conditions and thus do not contribute to the total energy or the equilibrium 
equations. We keep the first order terms in the model (2.9) since they do not sum to zero when the 
atomistic model is coupled to the continuum approximation in the quasicontinuum energy. The 



h 



i i jji 

+ 20F 



Uj — Uj-l 
h 



1 2 



Uj 



- u 

T 



j-2 



+ 



1 /// 



Uj 



- u 

T 



■3-2 



1 2 



(2.9) 
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atomistic energy (2.9) has the equilibrium equations 

(T a,h \ _ -<j>2F u 3+2 - <j>F u j+l + 2 (^F + ^2f) u 3 ~ $F u i~^ ~ <^2F u J-2 _ f 

[ U)j ~ h 2 ~ Jj: (2.10) 

Uj+2N = Uj, 

for — oo < j < oo. 

2.4. Linearized Continuum Model. For periodic u £ C°°(R) and Uj = u(xj), the total lin- 
earized atomistic energy 

S tot ' h (u) := S a ' h {u) - F(u) (2.11) 

converges to 

J [W(u'(x)) - f(x)u(x)} dx (2.12) 

as N — > oo, where the continuum strain energy density, W(e), is given by 

W(e) := [(4>' F + 20' 2F )e + + 40 2 ' F )e 2 ] • (2.13) 

The equilibrium equations (2.10) are a five-point consistent difference approximation of the equi- 
librium equation of the continuum model (2.12), which are 

-(# + 4&X = /. (2U) 
■u e (x + 2) = u e (x), 

for i£i 

Quasicontinuum approximations couple an approximation of the continuum model with the atom- 
istic model. The continuum approximation consists of a finite element discretization of the contin- 
uum model's elastic energy. The discretization uses a continuous, piecewise linear displacement u 
with the atom positions x as nodes. The external force term is applied as a point force at each 
node, so that (2.12) becomes 

N 

h[W(Dui)- fm]. (2.15) 

l=-N+l 

The continuum approximation has equilibrium equations 

- {(f) F + A<j) 2F ) — 2 =/,, -oo < I < oo, 

UI+2N = U h 

which is a three-point consistent difference approximation of the equilibrium equations for the 
continuum model (2.14). In one dimension, the above is actually the standard finite difference 
approximation of (2.12); however, it is framed in finite element terminology for flexibility in coars- 
ening, adapativity, and higher dimensional modelling. 

2.5. Assumptions. We assume that 

^ + 40' 2 V>O, (2.16) 

which implies that the total linearized atomistic energy (2.9) is positive definite (up to uniform 
translation of the displacement). Thus both equations (2.10) and (2.14) have a unique solution (up 
to uniform translation) provided that 

N 

£ /i = °- ( 2 - 17 ) 

j=-N+l 
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(2.18) 
(2.19) 

(2.20) 



For simplicity, we assume in the following that / is odd in addition to being periodic, that is, 

f(x) = —f(—x) and f(x + 2) = f(x) for — oo < x < oo, 
which implies that fj := f{xj) satisfies 

fj = -f-j and f j+2N = fj for - oo < j < oo. 
We obtain a unique, odd periodic solution satisfying the mean value condition 

N 

j=-N+l 

To give nonoscillatory solutions to the equilibrium equations (2.10) (that is, to guarantee that the 
roots of the corresponding characteristic equation are real (3.17)), we further assume that 

<f>p > and (j)2 F < 0. (2.21) 

The assumption (2.21) holds for potentials that allow an accurate second neighbor cut-off, such as 
the Lennard- Jones potential [9,10]. 

2.6. Energy-Based Quasicontinuum Approximation. The energy-based quasicontinuum ap- 
proximation (QCE) of £ a,h (u) decomposes the reference lattice into an atomistic region and a 
coarse-grained continuum region. It computes a total energy by using the atomistic energy (2.9) 
in the atomistic region and by using the continuum approximation (2.15) to sum the energy of the 
continuum region. 

For our analysis, we will consider an atomistic region defined by the atoms with reference 
positions Xj for j = —K,...,K, and a continuum region containing the remaining atoms, j = 
—N + 1, . . . , —K — 1 and j = K + 1, . . . , N. All atoms in the continuum region, along with the two 
atoms on the boundary, j = ±K will act as nodes for the continuum approximation. The continuum 
region can be decomposed into elements (x(/_i), xi) for / = — N + 1, • • ■ — K and / = K + 1, . . . , N. 
(In general, elements can contain many atoms of the reference lattice, but in this paper we do not 
consider coarsening in the continuum region.) 

To construct the contribution of the atomistic region to the total quasicontinuum energy, it is 
convenient to construct an energy associated with each atom by splitting equally the energy of each 
bond to obtain 



(u) := 





~U j+1 -Uj~ 




~U j+1 -Uj~ 


4>'f 


h 


h 



+ 



^F 









2 " 


~U j+2 -Uj~ 




'Uj+2-Uj' 




h 


h 





+- 



u 3 ~ u j~l 
h 



i 1 aJi 
+ 20F 



u 3 ~ U J~l 
h 



(2.22) 



+ 



^F 



~Uj - Uj- 2 ~ 


1 x» 


~Uj -Uj- 2 ~ 


to 


h 


h 





The continuum energy (2.15) is split into energy per element hW(Dui) where W is given in (2.13), 
and h = x\ — x\_\ is the length of the continuum element (x/_i,x^). 

To construct a quasicontinuum approximation QCE that conserves exactly the energy of atom- 
istic model (2.9) for lattices yj given by a uniform deformation gradient F (see (2.1)) the elements 
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(x-k-i, X-k) and (xk, xk+i) on the border of the atomistic region should contribute only one 
half of the continuum energy associated with that element. The QCE energy is then 

-K-l K 

£ ice,h^ ._ hW(Dm) + \hW{Du- K ) + J2 £ t k ( u ) 

(2-23) 

+ \hW{Du K+ i)+ Yl hw ( Du i)- 

l=K+2 

The equilibrium equations for the total QCE energy, £ qce > h (u) — F(u), then take the form [9, 10] 



L qce,h 



qce 



(2.24) 



where, for < j < N, we have 



qce,h 



U 



J-j+1 



2iij — v,i— 



7-1 



h 2 



+ < 



f ^2F 


-Uj+2 + 2Uj 


-Uj. 


-2 


Ah 2 






*<t>2F 


-Uj+2 + 2Uj 


- Uj- 


-2 


Ah 2 






H2F 


~Uj + 2 + 2Uj 


-Uj. 


-2 


Ah 2 








— Uj+i + 2uj 


-Uj. 


-1 


h 2 








—Uj + \ + 2uj 


-Uj. 


-1 


h 2 








—Uj + i + 2uj 


-Uj. 


-1 


h 2 







+ 



+ 



02F U j+2 ~ Uj 








h 2h ' 








2</>2F U j+l - U j 




4>2F 


Uj+2 ~ Uj 


h h 


h 


2h 


2(f>2 F Uj — Uj-i 


+ 


4>2F 


Uj - Uj-2 


h h 


h 


2h 


(f>2 F Uj - Uj-2 








h 2h ' 









0<j<K-2, 
j = K-l, 
j = K, 
j = K + l, 
j = K + 2, 
K + 3<j<N, 



with g given by 



9j 



0, < j < K - 2, 

3=K-1, 
J = K, 
j = K + l, 

0, K + 3 < j < N. 



"hi^Fi 

~2H§2Fi 
2h < r2Fi 



(2.25) 



For space reasons, we only list the entries for < j < N. The equations for all other j £ Z 
follow from symmetry and periodicity. Due to the symmetry in the definition of the atomistic and 
continuum regions, we have that Lff' h = L q ^^- and gj = —g~j for — N + 1 < i, j < 0. To see this, 



we define the involution operator (Su)j 
follows from the chain rule that 



-h-3 

-U- 



and observe that £ qce ' h (Su) = £i ce ' h (u). It then 



L qce ' h u — g - f for all periodic u. 



S T L gce,h Su _ s T g _ s T f 

Since S T = S and the assumption (2.19) is equivalent to 5/ = /, we can conclude that 

SL qce,h s = L qce,h wd S g = 9 . (2.26) 

Furthermore, we can conclude that the unique mean zero solution (2.20) to the equilibrium equa- 
tions (2.24) is odd. This follows from S -1 = S and (2.26) which together imply that Su is a 
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solution if and only if u is. Because S preserves the mean zero property, we conclude that u qce is 
odd. 

2.7. Stability of the Quasicontinuum Operator. Our analysis of the QCE error will utilize 
the following stability results for the operator L qce ' h . 



Lemma 2.1. If u := <\>" F — 5|</> 2 ' F | > 0, then 

hv-L qce > h v > v\\Dv\\% . 

h 

Proof. The stability result (2.27) follows from the identity 

-K-l K 

\hv ■ L qce K = hW(Dvi) + \hW{Dv- K ) + Y £j' H ( v ) 



(2.27) 



l=-N+l 



j=-K 



A 



+ \hW{Dv K+ i) + Y hW ( Dv i)i 

l=K+2 



where 



f«' h ( V ) 



h 




'V j+1 -Vj~ 


+ \$2F 


'vj+2-Vj' 


2 ' 


2 




h 




h 





+ 



1 Jjl 


'Vj-Vj-l' 


2 

i 1 x» 
+ 2V2F 


'vj-Vj- 2 ~ 


2 ' 


2<PF 


h 


h 





(2.28) 



and 



W(e) : =I(^ + 4^ F )6 2 . 



(2.29) 



We then have that 



N 



-K-l 



hvL^ h v>h Y h<t>F[(Dv j+1 ) 2 + {D Vj ) 2 ]-h Y 4 \<t>2F\(Dvj) 2 

j=-N+l j=-N+l 

K 

- 2h\$ F \(Dv- K ) 2 -hY \<%f\ [(Dv J+2 ) 2 + {Dv 3+1 ) 2 + {Dv 3 ) 2 + (Dvj^) 

j=-K 
N 

- 2h\<p'* F \(Dv K+1 ) 2 - h Y M^fKDvj) 2 

j=K+2 



>{4>"F-^\<t>2F 



A 



h y ( Dv rf 

j=-N+l 



□ 



The preceding stability Lemma 2.1 and the discrete Poincare inequality, 



|v|L 2 < 



h 



1 N 
■ n^ v ll/2 < H-DvlLa if V v j = (2.30) 

Zsm 2 * j=-N+l 

for < h < 1, give the following stability result in the ||-L2 norm. The proof of (2.30) follows from 

h 

verifying that (2sin^)//i is the smallest eigenvalue of D T D. 
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Lemma 2.2. If v := <p F — 5\<p2 F \ an d 



L qce,h v = b 



where Ylf=-N+i bj = 0, then 



\Dv\\f2 < ^- ||b||^ 2 



1 

2^ 



(2.31) 
(2.32) 



Proof. The result (2.32) follows from taking the inner product of (2.31) with v and then using the 
positive definiteness inequality (2.27) and the Poincare inequality (2.30). □ 

2.8. Quasi-nonlocal Quasicontinuum Approximation. The quasi-nonlocal quasicontinuum 
approximation (QNL) is similar to the QCE approximation, but it modifies the interactions around 
the interface in order to remove g from the elastic force. The quasi-nonlocal atoms ±K, ±(K + 1) 
interact directly with any atoms in the atomistic region within the next nearest neighbor cut-off, but 
interact as in the continuum region with other all other atoms. That is, unlike the atomistic model 
and continuum approximation, the form of energy contributions for quasi-nonlocal atoms depends 
on the type (atomstic, continuum, or quasi-nonlocal) of the neighboring atoms. For example, the 
energy contribution for j = K is 



(<p' F + 205 



2F) 



UK+l ~ U K 
h 



+ 1(^+4^) 



UK+l ~ U K 
h 



1 2 



h 

+ 2 



(p' F 


UK ~ UK-l 


1 aJ> 


UK ~ UK-l 




h 


h 



+ ■ 



>2F 



UK ~ U K -2 


1 1 jjl 
+ 2.V2F 


UK ~ U K -2 


2 ' 


h 


h 





and the energy contribution for j = K + 1 is 



{4>' F + 24>' 2F ) 



UK+2 — UK+l 



h 

+ 2 



UK+l ~ U K 



h 



h 



+ 



+ W'F+^2F) 



UK+2 - UK +l 
~h~ 



1 2 



I A" 
2^F 



UK+l ~ U K 



h 



+ ■ 



>2F 



UK+l - UK-l 

h 



i 1 aJ> 



UK+l - UK-l 

h 



The QNL energy is then 

-K-2 -K K-l 

S^\ U ):= hW(D Ul ) + ±hW(Du_ K -i)+ Y £]' h (u) + Y Sf (u) 

l=-N+l j=-K-l ' j=-K+l 

K+l N 

+ E £ f ( u ) + \hW{Du K+2 ) + Y hW(D Ul ). 

j=K ' l=K+3 

The QNL equilibrium equations are 

L« nl > h u qnl = f , 



(2.33) 
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where 



3+ A J — . 0<j<K-l, 

, // - ; u . !-'■!■ ill ~ u 3+2 + 2tt J+l ~ u j ■ _ ts 

4( ?2F <P2F ^2 ' 3 ~ K i 



Ah 2 






-u j+2 + 2Uj 


-Uj 


-2 


Ah 2 






— Uj+i + 2uj 


-Uj 


-1 


h 2 






—Uj+i + 2uj 


-Uj 


-1 



A^ F - Uj+1 ^^ j - Uj -\ l\ - J .V. 

We note that the QNL energy satisfies the symmetry condition S qnl,h (Su) = £ qnl,h (u), so the QNL 
operator L qnl,h is defined for j < by the identity SL qnl,h S = L qnl,h . While we have successfully 
removed the ghost force terms g, QNL is also not a consistent approximation of the continuum 
equations (2.14) at the interfacial atoms, such as j = K and j = K + 1 above. We will give a more 
detailed analysis of the approximation at the interface in Section 4. 

Our analysis of the QNL error will utilize the following stability result for the operator L qnl,h . 

Lemma 2.3. // 1 < p < oo, v := 4>" F - ^\4>2f\ > °' and 

L qnl,h v = b 



where J2f=-N+i = 0i then 



h\-L qnl ' h \ > v\\Dv\\}2 



h 



1 „ „ (2 - 34) 

\\Dv \\p2 < — b \\p2 . 

Proof. The proof of the stability result (2.34) follows the proof of the stability results for the QCE 
approximation in Lemmas 2.1 and 2.2 with the appropriate modification. □ 

Remark 2.1. The basic formulation of the QNL method removes the ghost force terms only for 
second-neighbor interactions in the ID case. A longer-range matching method is proposed in [12] 
that removes ghost forces for longer-range interactions by extending the region near the interface 
which have special energies. 

In 2D and 3D, there are similar restrictions on the interaction length that QNL corrects. Ad- 
ditional ghost forces arise when the quasicontinuum energy is extended to allow coarsening in the 
continuum region [12]. 

3. CONVERGENCE OF THE ENERGY-BASED QUASICONTINUUM SOLUTION 

We now analyze the quasicontinuum error and obtain estimates for its convergence rate by 
splitting the residual into two parts. One portion contains the low order terms, has support only 
near the atomistic to continuum interface, and is oscillatory. The remainder is higher order, and 
its influence will be bounded using the stability results. We recall that the QCE solution, u qce , is 
an odd, periodic solution of 

i w V = g + f, (3.1) 
and the continuum model solution is an odd, periodic function u e (x) satisfying 

-(4> F + A^ F )u': = f, (3.2) 
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Let u e denote the vector satisfying Uj = u e (xj). We will now derive estimates for the quasicontinuum 
error e = u e - u qce . 

It follows from the QCE equilibrium equation (3.1) that 

L qce,h e = L qce,h Ue _ JJl^h^^ = JJl^h^ _ g _ f . 

We split the residual L qce > h e as 

L qce,h e ._ p + (T 

where p contains the three lowest-order residual error terms in the interface, 

'0, 0<j<K-2, 

T '21-' ' r '21'" h -I/2) \ ~ h ( t ) 2F U K+l/2> +2^^2F M j^+l/2^' 3 = ^ ~ ^' 



(3.3) 
(3.4) 



(jl4>2F + ( I ) 2F U 'k+1/2) Ti + h$2F u 'k+l/2i +M^ ) 2F U K'+l/2^' 3 ~ ^, 
~ (jz<t>2F + &2F U 'k+1/2) Ti ~ h^2F U K+l/2' + 2^ < ^2F U K'+l/2^' 3 = K + 1) 

( 



(3.5) 



\$2F 



I) 2F U 'k+1/2 



) \ + w 



if+1/2' +24^2F M j^+l/2^' J — + 2, 

if + 3 < j < iV, 



and /9j = —p-j. Although pj = 0(1/ h) in the interface j = K — 1, . . . , K + 2, we will prove that the 
effect of p on the error is small away from the interface because it oscillates and the lowest order 
terms cancel in the sum 

K+2 

Ap:= PJ = ^2F< + i/ 2 - (3-6) 

j=K-l 

The residual term p represents the inconsistency of the operator L qce,h as a second-order finite 
difference approximation of the differential equation (3.2). This inconsistency is located only in the 
interface because the models themselves are second-order approximations. 

The residual term <x accounts for the error in approximating the continuum model (2.14) by a 
second-order finite difference approximation. We can estimate the residual a from Taylor's Theorem 
to obtain 



\a\\ t 2 < Ch 2 



u 



(4) 



L 2 



(3.7) 



Note that since u e , f, and g are odd, and p was constructed to be odd, then a is odd as well. 
Therefore, we can split the error e as 

G Gp ~\- Gfy 

such that 



L qce ' h e f 
L qce ' h e n 



P, 



(3.8) 



3.1. Global Discretization Error, e^. We now have by the stability (2.27) of L qce,h and the 
estimate of the residual (3.7) that 



l-De^H/a < Ch 2 

h 



,( 4 ) 



L 2 



(3.9) 



We can extend the bound to 
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Lemma 3.1. For defined in (3.8), we have 



l e <r||^ 



< V^WDe^Wp < Ch 2 



\De r 



< < 



Ch 2 



,( 4 ) 



3,1 

Ch 2 + p 



L 2 
(4) 



L 2 



,(4) 

L 2 

1 < P < 2, 

2 < p < oo. 



Proof. We obtain the Poincare inequality [7] 

llvLoo < ll-DvLi < \/2||£)vL2 
for all odd periodic v from the identity 



(3.10) 
(3.11) 

(3.12) 



fSLiM^) if J > 0, 
Vj {-Y.Lj-iKDve) if J < 0, 

which gives the first inequality in (3.12). The second follows from Holder's inequality. We can then 
obtain the error estimate (3.10) for Heo-Loo from the Poincare inequality (3.12) and the bound (3.9). 

h 

The "inverse" estimate [7] 



IDvLoo < hr lj2 IIDvl 



(2 
<h 



for all periodic v, and the Holder estimates [26] 

WDv 



< 



2 ~p „ 
2 2 p \\D\\ 

|2/p 



I 11^2 || 



,1-2/p 



1 <P< 2, 

2 < p < oo, 



combine to prove (3.11) by taking v = e^. 



(3.13) 



(3.14) 



□ 



3.2. Interfacial coupling error, e p . In the following, we will bound the error, e p , by constructing 
and estimating an explicit odd solution of 

L qce ' h e p = p. (3.15) 

Since pj is zero for all j except j = ±{K — 1,K, K + 1,K + 2}, e p satisfies a second-order, 
homogeneous recurrence relation in the interior of the continuum region and a fourth-order, ho- 
mogeneous recurrence relation in the interior of the atomistic region. Therefore, e p j is linear for 
j > K + 3 or j < — K — 3, and it is the sum of a linear solution and exponential solution for 
— K + 2 < j < K — 2. The coefficients for these solutions are determined by the equations in the 
atomistic to continuum interface. 

The homogeneous atomistic difference scheme 

-<t>2F u j+2 ~ <t>F u j+l + {2<t>F + 2 4>2F) u j ~ 4>F u j-\ ~ <t>2F u j-2 = (3.16) 

has characteristic equation 

-^ F A 2 - <$A + (2^ + 2^ F ) - f F k~ l - ^ F A- 2 = 0, 

with roots 



1,1, A,-, 



(3.17) 



where 



2F 



-20" 



2F 
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Based on the assumptions on <p in (2.16)and (2.21) and we have that A > 1. We note that if (j) 2F were 
positive contrary to assumption (2.21), then A would be negative which would give an oscillatory 
error in the atomistic region. General solutions of the homogeneous atomistic equations (3.16) have 
the form Uj = C\ + C 2 hj + C^A- 7 + C^\~\ but seeking an odd solution reduces this to the form 



u ; 



C 2 hj + C 3 (\i - A"'). 
The odd solution of the approximate error equations is thus of the form 



mi/ij + j8^) , 0<j<K, 
m 2 hj - m 2 + e K +i, j = K + l, 
[m 2 hj-m 2 , K + 2<j<N, 



(3.18) 



where expressing the unknown e Pt K+i using a perturbation of the linear solution, <zk+i, simplifies 
the solution of the equilibrium equations. The four coefficients mi, m 2 , eic+i, and (3 can be found 
by satisfying the four equilibrium equations in the interface, j = K — 1, . . . , K + 2. Summing the 
equilibrium equations across the interface gives 



Ap 



K+2 K+2 

E Pi = E ( L<?ce *" e 

j=K-l j=K-l 

e p,K-l - Zp,K-2 
h 2 

(<t>" F + 402f! 



+ 



PJJ 



b 2F 



£p,K + e p ,K-i — e P ,K-2 — e p ,K-3 



4h 2 



£p,K+3 — £p,K+2 



h 2 



= ^F + H 2F ){—-—)- 



h h 

The cancellation of the exponential terms in the final equality holds because 
<PU\ K - + + 02f)(A^- 1 - A~^ +1 - \ K ~ 2 + \' K+2 ) + ^ F (-\ K - Z + A^+ 3 ) = 0, 

which can be seen by summing (3.16) with the homogeneous solution yj = —X J for j = —K + 
2, ... ,K — 2. Thus, we have from summing the equilibrium equations (3.15) across the interface 
that 



mi=m 2 + 



hAp 



<t>" F + H 2F 



(3.19) 



The equality (3.19) can be interpreted as saying that the interfacial residual p acts as a source 
/ = Ap in the continuum equations (3.2) at x = xk- 



Lemma 3.2. For e p defined in (3.8), we have that 



Il e p||^oo < Ch(l + \u' K+1 / 2 \ + /i|"u^- +1 /2l + h\ u K+i/2\)> 
\De p \\ el < Ch l /P(l + \u' K+1/2 \ + h\u" K+1/2 \ + h\u% +1/2 \), 



(3.20) 



where C > is independent of h, K, and p, 1 < p < 00. 
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Proof. We will set up the system of equations for the coefficients in (3.18) and bound the decay of 
the coefficients. We split the interface equations as (A + hB)~x. = b, where 



A = 










'2F 
'2F 
'2F 
} 2F 



'2F 



4>2FlK+l - \4>2 F lK-l 
4>FlK+l + 02F1K+2 + \<\>2FlK 
-<P" F 1K ~ 2<t>2FlK ~ \4>2FlK-\ 
-\$2FlK 



2F 

4f 



D 



2F 



(K + l)<t> F + (\K + 2)<% F 
-\K^ F 



-(K + i)4>" F 

Kcj> F 



\K + 3)^ F 



(§#-§ 



2F 



2F 



(3.21) 



x 



mi 
rri2 
e-K+i 
P 



b = h 2 



Pk-i 

PK 

PK+1 

PK+2 



Using the equality (3.19), we rewrite the above as (Ak + hB)x = b where 

<t>2FlK+l ~ \<f> 2F lK-\ 



A K 



B 



2 < f>2F 
5 a" 



■if/)" 
2 < V2F 



^ 2F + -^ FlK -2^ FlK 



h 


%F 







" 




m 2 








x = 


e-K+i 













\$2FlK-\ 
~\$2FlK 



(3.22) 



b = h 2 



PK-l 
PK+1 
PK+2 



Ap 



<>F + 4<t>2F 



-K4>" F -{\K-\)^ F 
-\K^ F 



and 7j = v 3 . We have omitted the second equation, as the full system is linearly dependent 
after the elimination of m\ by (3.19). 

We note that Ak, B, and b do not depend on h directly, though Ak may have indirect dependence 
if K scales with h. Therefore, we can neglect B for sufficiently small h provided that A~^~ exists 
and is bounded uniformly in K. The following lemma, proven in [11], gives such a bound for Ak- 

Lemma 3.3. For all K satisfying 2 < K < N — 2, the matrix Ak is nonsingular and ll^^- 1 ]! < C 
where C > is independent of K and h. 

Due to the definition of p (3.5) and Ap (3.6), we have that 

||b||*» < C(h + h\u' K+1/2 \ + h 2 \u" K+l/2 \ + /i 2 K +1/2 |). 



The \n'" 
ine \u> K+1 /2 



contribution from Ap does not have h 3 as coefficient since K may scale linearly with 
N = 1/h. In general, we only have that hK < 1. 

Applying Lemma 3.3, we see that x is 0(h), and by (3.19), so is x. From (3.18), we finally 
conclude (3.20). □ 
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3.3. Total error. Combining the estimates (3.11) and (3.10) given in Lemma 3.1 for e CT with the 
estimate (3.20) in Lemma 3.2 for e p , we obtain from the triangle inequality that 

Theorem 3.1. Let e denote the QCE error. Then for 1 < p < oo, 2 < K < N — 2, and h 

sufficiently small, the error can be bounded by 



\\e\lqo <Ch(l + \u' K+1/2 \ + h\u" K+1/2 \ + h\u% +1/2 \ + h 
\De\lgP < Ch 1/p (l + \u K+1/2 \ + h\u" K+l/2 \ + h\u% +1/2 \ + h 



(4) 



(4) 



)■ 



L 2 



We note that the above argument giving optimal order estimates in Theorem 3.1 only utilized 
the estimate Ap =0(1), rather than the optimal estimate Ap =0(h) given in (3.6). 

Although our theorems give optimal order rates of convergence, our assumptions on the required 
regularity on u e is not optimal. We have assumed for simplicity of exposition that 
and used the estimate (3.7). Lower order estimates for a such as 



(4) 



L 2 



< OO 



| cr 1 1 ^ 2 < Ch 



2-8 



,(4-) 



L 2 



2 < s < 4, 



can be used to reduce the regularity assumptions on u e and still obtain optimal rates of convergence. 
Assuming the full regularity on u e also makes possible the precise identification and removal of the 
lower order terms in the error by the modification of the atomistic-to-continuum coupling scheme. 

4. Convergence of the Quasi-nonlocal Quasicontinuum Solution 
For the quasi-nonlocal approximation, we split the residual as 



L qnl,h Ue _ f = p + (T 



(4.1) 



where 



P= S 



o, 

<P2F u K+l/2 2 ( f } 2F a K+l/2 

fk'l 7 ," _ Id," II 1 " 

( f J 2F u K+l/2 2^2F a K+l/2 



X $2F<^nK j = K+l 



0, 



< j < K- 1, 
/;. j = K, 
j = K 

K + 2 < j < N, 



(4.2) 



and where 



< Ch 2 



u 



(4) 



LP 



The residual maximum norm ||pLoo here is 0(1) as opposed to the energy-based quasicontinuum 

h 

which has a 0(1 /h) residual maximum norm. However, the sum of p is similarly 0(h), that is, 

Ap = -h^ F u% +1/2 . (4.3) 
A similar argument as in the QCE case follows. We split the error as 

where 



yqnly 
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The same arguments apply to give the bounds (3.11) and (3.10) on e a . Thus, we need to work 
through the modified argument to bound e p . Since p is non-zero only at j = ±{K, K + 1}, the odd 



solution e p has the form 



mihj + /?(— ^7T 



< j < K, 
K + 1 < j < N. 



\m 2 hj - m 2 , 
Summing across the interface again gives 

K+2 K+2 

a p : = E pi = E ( Lqce,he J 



(4.4) 



j=K-l 



j=K-l 



Thus, we have again that 



m\ = ni2 + 



hAp 



A 



B 



^+H' 2 y (4,5) 

We focus on the equations at j = K — 1,K, and K + 1 and split the interface equations as 
(A + hB)x = b, where 

4>2F (f>2FlK+l 
(f)p + 2(fi% F <p' F ^K+l + 4>2 F lK+2 + i>2FlK 
-<j)p - 3<f% F -<j) F 7K ~ WiplK - i>2FlK-\ 

(K + 1)<P>> F -(K+1)<% F 

(K + l)(d> F + <f>'> F ) -(K+l)(<P'p + ^p) 

-K((j>'p + 3^'p) + 4>2F K^'p + 3<% F ) - (f>2F 








(4.6) 



X 



mi 
m 2 

(3 



h 2 



PK-l 

PK 

PK+1 



Using the equality (4.5), we rewrite the above as Ak~x. = b where 

4>2F 4>2FlK+l 
tf'p + 2<P'> F <t>'p lK+ l + <P'i F lK+2 + 4>2FlK 



A 



K 



X 



m 2 




h = h 2 



PK-l 
PK 



+ h 2 



Ap 



+ 



'2F 



{K + lWip 
(K+l)(<j> F + ^ F ) 



(4.7) 



We have omitted the second equation, as the full system is linearly dependent after substitution of 
mi. We have that Ak has full rank and ||b||^°° < Ch 2 (\u'^ +1 ^ 2 \ + |^ +1 y 2 |), so that we obtain the 
following error estimate for the quasi-nonlocal approximation. 

Theorem 4.1. Let e be the solution to the quasi-nonlocal error equation (4.1). Then for 1 < p < oo, 
2 < K < N — 2, and h sufficiently small, the error can be bounded by 



u 



\\e\\ ( oo < Ch 2 (\u" K+1 i 2 \ + \u"k +1 j 2 \ + 

\\De\\ eh < Ch 1+1 ^ {\u" K+l/2 \ + |< +1/2 | + 
where C > is independent of h, K, and p. 



(4) 



L 2 



)• 
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We note that the proof above of the optimal order estimates for the quasi-nonlocal approximation 
does use the full 0(h) order of the estimate (4.3) for Ap. 
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